Removed DESeq part and added some file saving for CSVs
# Load libraries
library(DESeq2)
library(edgeR)
library(limma)
library(variancePartition)
library(foreach)
# library(clusterProfiler)
library(org.Hs.eg.db)
# library(ggh4x)
library(magrittr)
library(enrichplot)
library(patchwork)
library(fgsea)
library(doParallel)
library(tidyverse)
# Define main path
main_path <- file.path("/home/jbrenton/ASAP_Hardy_post_align_files/Final_Dataset/Hardy_Limma_Analysis_Final_version_only/")
# main_path <- here::here()
# Load helper functions
source(file.path(main_path, "R/hf_graph_and_themes.R"))
source(file.path(main_path, "R/hf_reports.R"))
# Set options
options(dplyr.summarise.inform = FALSE)
options(lifecycle_verbosity = "warning")
knitr::opts_chunk$set(echo = F, warning = T, message = T,
out.width="100%", fig.align = "center", dpi = 150,
fig.width = 7.2, fig.height = 7.2)## Custom sciRmdTheme with increased width.
## Please ask guillermorocamora@gmail.com for more details.
## Download from devtools::install_github("https://github.com/guillermo1996/grpSciRmdTheme")
if(require(grpSciRmdTheme)){
grpSciRmdTheme::set.theme(
theme = "default",
color = NULL,
header.sticky = FALSE,
list.group.icon = "arrow",
font.family = "Arial",
font.color = "black",
header.color = "darkblue",
figure.perc = 100
)
}## Custom styles for the output html
cat('
<style type="text/css">
.dataTables_scrollHeadInner{
width:100% !important;
}
.dataTables_scrollHeadInner table{
width:100% !important;
}
/*
.code-folding-btn {
display: none;
}
*/
h3, .h3 {
font-size: 22px!important;
}
h4, .h4 {
font-size: 18px!important;
}
h5, .h5 {
font-size: 16px!important;
}
body{
font-size: 13px;
}
.tocify-subheader {
text-indent: 15px;
display: none;
font-size: 12px;
}
.break5 .nav-tabs, .break4 .nav-tabs, .break3 .nav-tabs, .break2 .nav-tabs{
display: grid;
text-align: center;
font-size: 14px;
font-weight: bold;
}
.break5 .nav-tabs:before, .break4 .nav-tabs:before, .break3 .nav-tabs:before, .break2 .nav-tabs:before,
.break5 .nav-tabs:after, .break4 .nav-tabs:after, .break3 .nav-tabs:after, .break2 .nav-tabs:after {
content: unset;
}
.break5 .nav-tabs{
grid-template-columns: repeat(5, 1fr);
}
.break4 .nav-tabs{
grid-template-columns: repeat(4, 1fr);
}
.break3 .nav-tabs{
grid-template-columns: repeat(3, 1fr);
}
.break2 .nav-tabs{
grid-template-columns: repeat(2, 1fr);
}
.breakB5 .nav-tabs{
grid-template-columns: repeat(5, 1fr);
font-size: 12px;
}
.breakB4 .nav-tabs{
grid-template-columns: repeat(4, 1fr);
font-size: 12px;
}
.breakB3 .nav-tabs{
grid-template-columns: repeat(3, 1fr);
font-size: 12px;
}
.breakB2 .nav-tabs{
grid-template-columns: repeat(2, 1fr);
font-size: 12px;
}
.stop-break .nav-tabs {
display: flex;
font-size: 13px;
font-weight: normal;
}
.stop-break .nav-tabs:before,
.stop-break .nav-tabs:after {
content: unset;
}
</style>')This report includes a differential gene expression analysis on the Hardy bulk dataset with limma package, using voom and duplicateCorrelation. It is created as a template and not as a report with the final results.
The main path to run this report on the RytenLab server is at:
/home/grocamora/RytenLab-Research/18-DGE_limma_dream/
This report is meant to be run after after:
Scripts/01-Hardy_covariate_identification.RRMarkdown/02-Hardy_covariate_identification.RRMarkdown/03-PC_Analysis_Outliers.RmdThe most relevant data to load prior to the analysis are the metadata, the salmon quantifications, a list of blacklisted genes and the outliers to remove:
# Input paths
main_path <- "/home/jbrenton/ASAP_Hardy_post_align_files/Final_Dataset/Hardy_Limma_Analysis_Final_version_only/"
results_path <- file.path(main_path, "results")
txi.salmon_path <- file.path(main_path, "data/Hardy/txi.salmon.rds")
sample_metadata_path <- file.path(main_path, "metadata/Hardy/updated_metadata_w_PMI.csv")
metadata_path <- file.path(results_path, "all_brain_areas_covariates.rds")
sample_outliers_path <- file.path(results_path, "sample_outliers16.rds")
blacklist_genes <- file.path(main_path, "data/blist_genes_no_vers.csv")
gene_map_path <- file.path(main_path, "data/gencode_txid_to_geneid.txt")
# Output paths
hardy_dge_results_path <- file.path(results_path, "dge_variables/")
deg_counts_output_path <- file.path(hardy_dge_results_path, "PD_Collapsed_dge_counts.rds")
dge_all_path <- file.path(hardy_dge_results_path, "PD_Collapsed_dge_all.rds")
dir.create(hardy_dge_results_path, showWarnings = F, recursive = F)
# Load data
txi.salmon <- readRDS(txi.salmon_path)
cts <- txi.salmon$counts
metadata <- loadMetadataHardy(metadata_path, sample_metadata_path)
sample_outliers <- readRDS(sample_outliers_path)
bxp_outliers <- sample_outliers$bxp_id_full[sample_outliers$outlier]
# Load blacklisted genes
blacklist_genes <- vroom::vroom(blacklist_genes, col_names = T, delim = ",")
blacklist_genes <- as.vector(unlist(blacklist_genes))
gencode_txid_to_geneid <- vroom::vroom(gene_map_path) %>%
`colnames<-`(c("tx_id", "gene_id","gene_name", "description")) %>%
dplyr::mutate(tx_id = sub("\\..+", "", tx_id),
gene_id = sub("\\..+", "", gene_id))
# Recode PD to combine PD groups
metadata %<>% mutate(pd=case_when(pd=="PD" ~ "PD",
pd=="PDD" ~ "PD",
pd=="control" ~ "control"))The covariates employed to model the gene expression data are extracted from the Covariate identification report previously mentioned and are the following:
GC_NC_40_59, PCT_UTR_BASES, gender, Total_Sequences, INTRONIC_BASES and tss_up_1kb_tag_pct.
The experimental condition of interest in this report will be the type of PD of the samples: controls, PD without dementia (PD) and PD with dementia (PDD). We will also study two different scenarios for the brain regions: i) a collapsed brain region approach and ii) a per brain region analysis.
It is also required to fill in certain parameters that will affect the pipeline and the results:
Updated to final decision on outliers and gene expression threshold.
Removed Gene Enrichment Analyses from this report.
Added GO, KEGG and REACTOME gene enrichment analyses.
Updated default report to include 16 outliers with the DGEreport::degFilter threshold.
Merged the two limma pipelines to combine diagnosis and brain region. Common practice in the literature.
Added DEGreport filter to the gene expression threhsolds.
Added a summary table of all results.
Removed the dream pipeline.
Added the limma/voom pipeline for all tissues combined and by brain region.
Added DESeq2 analysis for reference.
Added additional study about the gene expression threshold.
limma/voom and duplicateCorrelationA complete tutorial on how to use limma can be found in their documentation page. On its core, limma fits a linear model to assess differential gene expression across different comparisons (or contrasts).
Given that we are comparing across different tissues with samples from the same individuals, we need to treat the sample donor as a random effect. In limma, this is handled by the duplicateCorrelation function, which estimates a correlation between replicates from the same individual and gene expressions. In this approach, a single correlation value is employed genome-wide, assuming that the correlation is the same across all genes. Other alternatives, such as the dream package, estimate a per-gene correlation, but their study is beyond the scope of this report.
Following the instructions in the official documentation and in this tutorial about design matrices, the two factors of interests (PD status and brain region) are grouped into a single factor. Thus, only one model is fitted and all comparisons are extracted from it.
The first step is to generate a dataframe with the covariates of interests (limma_covs) and a dataframe with the covariates plus the grouping factors (limma_covs_all) which includes the experimental condition, the brain region and the individual ID. Outliers are also removed following the criteria explained in PCA Analysis for outliers report: outliers detected using PCA Z-scores and WGCNA sample connectivity Z-scores. It is important to ensure that the order of samples in the counts object is the same as the order of samples in the metadata table.
# Ensure sample order
metadata <- metadata[sapply(metadata$bxp_id_full, function(x) grep(x, colnames(cts))), ]
stopifnot(all(colnames(cts) == metadata$bxp_id_full))
# Covariate + grouping dataframe
limma_covs_all <- metadata %>%
dplyr::select(sample_id, bxp_id_full, individual_id = case_id,
condition = all_of(experimental_condition),
region = all_of(experimental_region),
all_of(covariates)) %>%
# Filter the outliers
dplyr::filter(!bxp_id_full %in% bxp_outliers) %>%
# Generate the grouping factor as "condition_region"
dplyr::mutate(group = as.factor(paste0(condition, "_", region))) %>%
dplyr::relocate(individual_id, group) %>%
# Convert all characters to factor
dplyr::mutate(across(where(is.character), ~ as.factor(.))) %>%
# Scale the covariates again after the removal of outliers
dplyr::mutate(across(where(is.numeric), ~ scale(., center = T, scale = T))) %>%
tibble::column_to_rownames("bxp_id_full")
# Covariates dataframe
limma_covs <- limma_covs_all %>% dplyr::select(all_of(covariates))Next, we generate the DGE object and apply the following preprocessing steps:
DEGreport: employs the function DEGreport::degFilter, which requires at least one count in all samples of at least one group.edgeR: employs the function edgeR::filterByExpr which, in this situation, requires a minimum cpm of 0.11 in at least 19 samples.cpm: requires a minimum cpm of one in at least 50% of the samples.# The author recommends to use "keep.lib.sizes = FALSE" after gene filtering,
# but before "calcNormFactors". Source:
# https://support.bioconductor.org/p/9136787/
# Remove outliers
cts_clean <- cts[, !colnames(cts) %in% bxp_outliers]
# Generate the DGE object
dge0 <- edgeR::DGEList(cts_clean, group = limma_covs_all$group)
# Remove the blacklisted genes
dge1 <- dge0[!rownames(dge0) %in% blacklist_genes, , keep.lib.sizes = F]
# Remove low-expressed genes. We are currently testing different approaches to
# filtering the low-expressed genes. Three valid approaches are currently
# employed:
#
# i) edgeR: use the function "edgeR::filterByExpr". It currently requires a cpm
# > 0.11 in at least 19 samples.
#
# ii) cpm: the same filter applied to covariate identification and PC analysis.
# It requires a cpm or at least 1 in 50% of the samples.
#
# iii) DEGreport: a previously used filter by the group. It requires at least
# one count across all samples in one of the groups.
if(gene_filter_method == "DEGreport"){
filter_count <- DEGreport::degFilter(counts = dge1$counts,
metadata = limma_covs_all,
group = "group",
min = 1, # All samples in group must have more than expr > 0
minreads = 0)
isexpr <- rownames(dge1$counts) %in% rownames(filter_count)
}else if(gene_filter_method == "edgeR"){
isexpr <- edgeR::filterByExpr(dge1, group = limma_covs_all$group)
}else if(gene_filter_method == "cpm"){
isexpr <- rowSums(edgeR::cpm(dge1) > 1) >= 0.5*ncol(dge1)
}else{
isexpr <- rep(T, nrow(dge1$counts))
}
dge2 <- dge1[isexpr, , keep.lib.sizes = F]
# Calculate normalization factors - adds that column to the samples list
dge3 <- edgeR::calcNormFactors(dge2)
# Final DGE object
dge <- dge3
# Clear variables
rm(cts_clean, dge0, dge1, dge2, dge3)
message(paste0("Employed gene expression filter: ", gene_filter_method)){bg.success} ## Employed gene expression filter: DEGreport
Finally, we generate the design object required for the limma pipeline. We also specify the formula that the linear model fit will use, and it must be shaped as ~ 0 + experimental_variable + covariates. The use of intercept is advised to allow for all different comparisons. The formula employed in the analysis is the following:
Formula: ~ 0 + group + gender + GC_NC_40_59 + tss_up_1kb_tag_pct + PCT_UTR_BASES + Total_Sequences + INTRONIC_BASES
limma_formula <- reformulate(c(0, "group", limma_covs %>% dplyr::select(where(is.factor), where(is.numeric)) %>% colnames()))
limma_factor_formula <- reformulate(c(0, "group", limma_covs %>% dplyr::select(where(is.factor)) %>% colnames()))
limma_num_formula <- reformulate(limma_covs %>% dplyr::select(where(is.numeric)) %>% colnames())
limma_mod <- model.matrix(limma_factor_formula, data = limma_covs_all) %>%
`colnames<-`(stringr::str_remove_all(colnames(.), "^group"))
limma_design <- data.frame(limma_mod, limma_covs %>% dplyr::select(where(is.numeric)))duplicateCorrelationWe use voom to obtain the weights of each gene based on the residuals of a linear fit (more information). Since we also want to consider the correlation between samples from the same individual, we employ duplicateCorrelation. Both functions are applied twice based on recommendations from the developers.
By default, the following code block stores and load from memory the vobj and dupcor variables. Please remove files from disk or modify the results folder to run the analysis again.
vobj_path <- file.path(hardy_dge_results_path, "PD_Collapsed_vobj_limma.rds")
dupcor_path <- file.path(hardy_dge_results_path, "PD_Collapsed_dupcor_limma.rds")
if(!file_test("-f", vobj_path, dupcor_path)){
vobj_tmp <- limma::voom(dge, limma_design, save.plot = T)
dupcor_tmp <- limma::duplicateCorrelation(vobj_tmp, limma_design, block = limma_covs_all$individual_id)
vobj <- limma::voom(dge, limma_design, save.plot = T,
block = limma_covs_all$individual_id, correlation = dupcor_tmp$consensus)
dupcor <- limma::duplicateCorrelation(vobj, limma_design, block = limma_covs_all$individual_id)
vobj %>% saveRDS(vobj_path)
dupcor %>% saveRDS(dupcor_path)
}else{
vobj <- readRDS(vobj_path)
dupcor <- readRDS(dupcor_path)
}
plotVobj(vobj)The next step is to fit the model for each gene. We employ the function lmFit followed by eBayes to estimate the log-odds of differential expression and t-statistic.
fit_path <- file.path(hardy_dge_results_path, "PD_Collapsed_fit_limma.rds")
if(!file_test("-f", fit_path)){
fit <- limma::lmFit(vobj, limma_design, block = limma_covs_all$individual_id, correlation = dupcor$consensus)
fit <- limma::eBayes(fit)
fit %>% saveRDS(fit_path)
}else{
fit <- readRDS(fit_path)
}One of the most convenient aspects of the limma pipeline is the ability to accommodate arbitrary experimental complexity. The function limma::contrasts.fit allows for comparisons between groups to be made after the fit. To generate these contrasts, we need to define a contrast matrix specifing the groups to compare. It is important to notice that each positive contribution must add to 1 and each negative contribution must add to -1 (otherwise, the log fold-changes are scaled). The reference group must have a negative contribution.
contrasts_fit <- limma::makeContrasts(
"Collapse-PD-Control" = (PD_ACG + PD_IPL + PD_MFG + PD_MTG -
control_ACG - control_IPL - control_MFG - control_MTG)/4,
"ACG-PD-Control" = PD_ACG - control_ACG,
"IPL-PD-Control" = PD_IPL - control_IPL,
"MFG-PD-Control" = PD_MFG - control_MFG,
"MTG-PD-Control" = PD_MTG - control_MTG,
levels = limma_design)
# Fit by the contrasts
fit_con <- limma::eBayes(limma::contrasts.fit(fit, contrasts_fit))## Warning in limma::contrasts.fit(fit, contrasts_fit): row names of contrasts
## don't match col names of coefficients
In our scenario, a total of 15 different contrasts were made: 3 different case vs. control studies (i.e. PD vs control, PDD vs control and PDD vs PD) in 5 different tissue combinations (collapsed and by brain region). The following table summarizes the contrasts:
In this section, we will compare the results of combining all brain regions. The generated contrast table is drawn below, where each coefficient is set to \(1/\text{number of brain regions}\). As mentioned previously, the sign of the coefficient is important, leaving the negative contributions to the reference group.
To extract the results from the fit object, we employ a custom function that internally calls limma::topTable for a given coefficient (i.e. the name for the contrasts of interest in the contrasts matrix):
dge_collapsed_i <- extractResultsDGE(fit_con, coef = "Collapse-PD-Control")
dge_collapsed <- dplyr::bind_rows(
dge_collapsed_i %>% dplyr::mutate(tissue = "Collapse", test = "PD-Control"))The number of differentially expressed genes obtained for each comparison is summarized in the following table:
## Down NotSig Up
## 123 22786 62
Volcano plots allow for a more global view of the DGE results by showing statistical significance vs magnitude of change - i.e. log transformed p-values in the Y-axis and the log2 fold-change in the X-axis. Our requirement for significance is \(p-value <= 0.05\), drawn as a red line. Every gene above the mentioned line is considered differentially expressed. The sign of the log2 fold-change specifies whether it is upregulated or downregulated. Additionally, the top most significant differentially expressed genes are named in each test.
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text_repel()`).
Here we represent the 10 most significant differentially expressed genes along with their magnitude of change and their statistical significance - i.e. log transformed p-values as the fill colour and log2 fold-change in the Y-axis.
The mean-difference plots are another useful tool to visualize the DGE results. Instead of showing the statistical significance (as the Volcano plot), they represent the magnitude of change against the average log-CPM expression - i.e. average log-CPM in the x-Axis and log2 fold-change in the X axis. Significant differentially expressed genes are represented both by their color and their shape:
# Extract the results from the main fit.
tissue_contrasts <- colnames(contrasts_fit)[-grep("Collapse", colnames(contrasts_fit))]
dge_tissues <- lapply(tissue_contrasts, function(contrast_name){
tissue <- stringr::str_split_fixed(contrast_name, "-", n = 2)[1]
test <- stringr::str_split_fixed(contrast_name, "-", n = 2)[2]
extractResultsDGE(fit_con, coef = contrast_name) %>% dplyr::mutate(tissue = tissue, test = test)
}) %>% dplyr::bind_rows()In this section, we will compare the results for the different brain regions.
## Warning: 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")
| Tissue - Test | Downregulated | Unchanged | Upregulated |
|---|---|---|---|
| ACG | |||
| PD-Control | 48 | 22909 | 14 |
| IPL | |||
| PD-Control | 17 | 22952 | 2 |
| MFG | |||
| PD-Control | 33 | 22932 | 6 |
| MTG | |||
| PD-Control | 49 | 22870 | 52 |
In this section, a brief summary of the results will be presented.
First, we can study the number of differentially expressed genes, regardless of their magnitude or direction of change:
Each tile contains the proportion of upregulated/downregulated genes. Only proportions higher than 30% are written in the tiles for clarity:
## Warning: 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")
| Tissue | Regulated | Number of DEG - Limma |
|---|---|---|
| Collapse | ||
| Collapse | Significantly Downregulated | 123 |
| Collapse | Significantly Upregulated | 62 |
| ACG | ||
| ACG | Significantly Downregulated | 48 |
| ACG | Significantly Upregulated | 14 |
| IPL | ||
| IPL | Significantly Downregulated | 17 |
| IPL | Significantly Upregulated | 2 |
| MFG | ||
| MFG | Significantly Downregulated | 33 |
| MFG | Significantly Upregulated | 6 |
| MTG | ||
| MTG | Significantly Downregulated | 49 |
| MTG | Significantly Upregulated | 52 |
## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.3.0 (2023-04-21)
## os Ubuntu 20.04.6 LTS
## system x86_64, linux-gnu
## ui X11
## language (EN)
## collate en_GB.UTF-8
## ctype en_GB.UTF-8
## tz Etc/UTC
## date 2025-09-30
## pandoc 3.1.1 @ /usr/lib/rstudio-server/bin/quarto/bin/tools/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## abind 1.4-8 2024-09-12 [1] CRAN (R 4.3.0)
## AnnotationDbi * 1.64.1 2023-11-03 [1] Bioconductor
## aod 1.3.3 2023-12-13 [1] CRAN (R 4.3.0)
## ape 5.8-1 2024-12-16 [1] CRAN (R 4.3.0)
## aplot 0.2.4 2024-12-17 [1] CRAN (R 4.3.0)
## backports 1.5.0 2024-05-23 [1] CRAN (R 4.3.0)
## Biobase * 2.62.0 2023-10-24 [1] Bioconductor
## BiocGenerics * 0.48.1 2023-11-01 [1] Bioconductor
## BiocParallel * 1.36.0 2023-10-24 [1] Bioconductor
## Biostrings 2.70.3 2024-03-13 [1] Bioconductor 3.18 (R 4.3.0)
## bit 4.6.0 2025-03-06 [1] CRAN (R 4.3.0)
## bit64 4.6.0-1 2025-01-16 [1] CRAN (R 4.3.0)
## bitops 1.0-9 2024-10-03 [1] CRAN (R 4.3.0)
## blob 1.2.4 2023-03-17 [1] CRAN (R 4.3.0)
## bookdown 0.42 2025-01-07 [1] CRAN (R 4.3.0)
## boot 1.3-31 2024-08-28 [1] CRAN (R 4.3.0)
## broom 1.0.7 2024-09-26 [1] CRAN (R 4.3.0)
## bslib 0.8.0 2024-07-29 [1] CRAN (R 4.3.0)
## cachem 1.1.0 2024-05-16 [1] CRAN (R 4.3.0)
## caTools 1.18.3 2024-09-04 [1] CRAN (R 4.3.0)
## circlize 0.4.16 2024-02-20 [1] CRAN (R 4.3.0)
## cli 3.6.5 2025-04-23 [1] CRAN (R 4.3.0)
## clue 0.3-66 2024-11-13 [1] CRAN (R 4.3.0)
## cluster 2.1.8.1 2025-03-12 [1] CRAN (R 4.3.0)
## codetools 0.2-20 2024-03-31 [1] CRAN (R 4.3.0)
## colorspace 2.1-1 2024-07-26 [1] CRAN (R 4.3.0)
## ComplexHeatmap 2.18.0 2023-10-24 [1] Bioconductor
## ConsensusClusterPlus 1.66.0 2023-10-24 [1] Bioconductor
## corpcor 1.6.10 2021-09-16 [1] CRAN (R 4.3.0)
## cowplot 1.1.3 2024-01-22 [1] CRAN (R 4.3.0)
## crayon 1.5.3 2024-06-20 [1] CRAN (R 4.3.0)
## crosstalk 1.2.1 2023-11-23 [1] CRAN (R 4.3.0)
## data.table 1.17.6 2025-06-17 [1] CRAN (R 4.3.0)
## DBI 1.2.3 2024-06-02 [1] CRAN (R 4.3.0)
## DEGreport 1.38.5 2023-12-06 [1] Bioconductor 3.18 (R 4.3.0)
## DelayedArray 0.28.0 2023-10-24 [1] Bioconductor
## DESeq2 * 1.42.1 2024-03-06 [1] Bioconductor 3.18 (R 4.3.0)
## dichromat 2.0-0.1 2022-05-02 [1] CRAN (R 4.3.0)
## digest 0.6.37 2024-08-19 [1] CRAN (R 4.3.0)
## doParallel * 1.0.17 2022-02-07 [1] CRAN (R 4.3.0)
## DOSE 3.28.2 2023-12-10 [1] Bioconductor 3.18 (R 4.3.0)
## dplyr * 1.1.4 2023-11-17 [1] CRAN (R 4.3.0)
## DT 0.33 2024-04-04 [1] CRAN (R 4.3.0)
## edgeR * 4.0.16 2024-02-18 [1] Bioconductor 3.18 (R 4.3.0)
## enrichplot * 1.22.0 2023-10-24 [1] Bioconductor
## EnvStats 3.0.0 2024-08-24 [1] CRAN (R 4.3.0)
## evaluate 1.0.4 2025-06-18 [1] CRAN (R 4.3.0)
## fANCOVA 0.6-1 2020-11-13 [1] CRAN (R 4.3.0)
## farver 2.1.2 2024-05-13 [1] CRAN (R 4.3.0)
## fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.3.0)
## fastmatch 1.1-6 2024-12-23 [1] CRAN (R 4.3.0)
## fgsea * 1.28.0 2023-10-24 [1] Bioconductor
## forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.3.0)
## foreach * 1.5.2 2022-02-02 [1] CRAN (R 4.3.0)
## fs 1.6.6 2025-04-12 [1] CRAN (R 4.3.0)
## generics 0.1.4 2025-05-09 [1] CRAN (R 4.3.0)
## GenomeInfoDb * 1.38.8 2024-03-15 [1] Bioconductor 3.18 (R 4.3.0)
## GenomeInfoDbData 1.2.11 2023-11-20 [1] Bioconductor
## GenomicRanges * 1.54.1 2023-10-29 [1] Bioconductor
## GetoptLong 1.0.5 2020-12-15 [1] CRAN (R 4.3.0)
## ggdendro 0.2.0 2024-02-23 [1] CRAN (R 4.3.0)
## ggforce 0.4.2 2024-02-19 [1] CRAN (R 4.3.0)
## ggfun 0.1.8 2024-12-03 [1] CRAN (R 4.3.0)
## ggplot2 * 4.0.0 2025-09-11 [1] CRAN (R 4.3.0)
## ggplotify 0.1.2 2023-08-09 [1] CRAN (R 4.3.0)
## ggraph 2.2.1 2024-03-07 [1] CRAN (R 4.3.0)
## ggrepel 0.9.6 2024-09-07 [1] CRAN (R 4.3.0)
## ggtree 3.10.1 2024-02-25 [1] Bioconductor 3.18 (R 4.3.0)
## GlobalOptions 0.1.2 2020-06-10 [1] CRAN (R 4.3.0)
## glue 1.8.0 2024-09-30 [1] CRAN (R 4.3.0)
## GO.db 3.18.0 2023-11-23 [1] Bioconductor
## GOSemSim 2.28.1 2024-01-17 [1] Bioconductor 3.18 (R 4.3.0)
## gplots 3.2.0 2024-10-05 [1] CRAN (R 4.3.0)
## graphlayouts 1.2.1 2024-11-18 [1] CRAN (R 4.3.0)
## gridExtra 2.3 2017-09-09 [1] CRAN (R 4.3.0)
## gridGraphics 0.5-1 2020-12-13 [1] CRAN (R 4.3.0)
## grpSciRmdTheme * 0.4.0.1 2023-12-03 [1] Github (guillermo1996/grpSciRmdTheme@5f7645e)
## gtable 0.3.6 2024-10-25 [1] CRAN (R 4.3.0)
## gtools 3.9.5 2023-11-20 [1] CRAN (R 4.3.0)
## HDO.db 0.99.1 2023-11-23 [1] Bioconductor
## hms 1.1.3 2023-03-21 [1] CRAN (R 4.3.0)
## htmltools 0.5.8.1 2024-04-04 [1] CRAN (R 4.3.0)
## htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.3.0)
## httr 1.4.7 2023-08-15 [1] CRAN (R 4.3.0)
## igraph 2.1.4 2025-01-23 [1] CRAN (R 4.3.0)
## IRanges * 2.36.0 2023-10-24 [1] Bioconductor
## iterators * 1.0.14 2022-02-05 [1] CRAN (R 4.3.0)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.3.0)
## jsonlite 2.0.0 2025-03-27 [1] CRAN (R 4.3.0)
## kableExtra 1.4.0 2024-01-24 [1] CRAN (R 4.3.0)
## KEGGREST 1.42.0 2023-10-24 [1] Bioconductor
## KernSmooth 2.23-26 2025-01-01 [1] CRAN (R 4.3.0)
## knitr 1.49 2024-11-08 [1] CRAN (R 4.3.0)
## labeling 0.4.3 2023-08-29 [1] CRAN (R 4.3.0)
## lattice 0.22-7 2025-04-02 [1] CRAN (R 4.3.0)
## lazyeval 0.2.2 2019-03-15 [1] CRAN (R 4.3.0)
## lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.3.0)
## limma * 3.58.1 2023-10-31 [1] Bioconductor
## lme4 1.1-36 2025-01-11 [1] CRAN (R 4.3.0)
## lmerTest 3.1-3 2020-10-23 [1] CRAN (R 4.3.0)
## locfit 1.5-9.10 2024-06-24 [1] CRAN (R 4.3.0)
## logging 0.10-108 2019-07-14 [1] CRAN (R 4.3.0)
## lubridate * 1.9.4 2024-12-08 [1] CRAN (R 4.3.0)
## magrittr * 2.0.4 2025-09-12 [1] CRAN (R 4.3.0)
## MASS 7.3-60.0.1 2024-01-13 [1] CRAN (R 4.3.0)
## Matrix 1.6-5 2024-01-11 [1] CRAN (R 4.3.0)
## MatrixGenerics * 1.14.0 2023-10-24 [1] Bioconductor
## matrixStats * 1.5.0 2025-01-07 [1] CRAN (R 4.3.0)
## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.3.0)
## minqa 1.2.8 2024-08-17 [1] CRAN (R 4.3.0)
## mnormt 2.1.1 2022-09-26 [1] CRAN (R 4.3.0)
## mvtnorm 1.3-3 2025-01-10 [1] CRAN (R 4.3.0)
## nlme 3.1-166 2024-08-14 [1] CRAN (R 4.3.0)
## nloptr 2.2.1 2025-03-17 [1] CRAN (R 4.3.0)
## numDeriv 2016.8-1.1 2019-06-06 [1] CRAN (R 4.3.0)
## org.Hs.eg.db * 3.18.0 2023-11-23 [1] Bioconductor
## patchwork * 1.3.0 2024-09-16 [1] CRAN (R 4.3.0)
## pbkrtest 0.5.3 2024-06-26 [1] CRAN (R 4.3.0)
## pillar 1.11.0 2025-07-04 [1] CRAN (R 4.3.0)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.0)
## plyr 1.8.9 2023-10-02 [1] CRAN (R 4.3.0)
## png 0.1-8 2022-11-29 [1] CRAN (R 4.3.0)
## polyclip 1.10-7 2024-07-23 [1] CRAN (R 4.3.0)
## psych 2.4.12 2024-12-23 [1] CRAN (R 4.3.0)
## purrr * 1.1.0 2025-07-10 [1] CRAN (R 4.3.0)
## qvalue 2.34.0 2023-10-24 [1] Bioconductor
## R6 2.6.1 2025-02-15 [1] CRAN (R 4.3.0)
## rbibutils 2.3 2024-10-04 [1] CRAN (R 4.3.0)
## RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.3.0)
## Rcpp 1.1.0 2025-07-02 [1] CRAN (R 4.3.0)
## RCurl 1.98-1.17 2025-03-22 [1] CRAN (R 4.3.0)
## Rdpack 2.6.4 2025-04-09 [1] CRAN (R 4.3.0)
## readr * 2.1.5 2024-01-10 [1] CRAN (R 4.3.0)
## reformulas 0.4.0 2024-11-03 [1] CRAN (R 4.3.0)
## remaCor 0.0.18 2024-02-08 [1] CRAN (R 4.3.0)
## reshape 0.8.10 2025-06-19 [1] CRAN (R 4.3.0)
## reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.3.0)
## RhpcBLASctl 0.23-42 2023-02-11 [1] CRAN (R 4.3.0)
## rjson 0.2.23 2024-09-16 [1] CRAN (R 4.3.0)
## rlang 1.1.6 2025-04-11 [1] CRAN (R 4.3.0)
## rmarkdown 2.29 2024-11-04 [1] CRAN (R 4.3.0)
## RSQLite 2.3.9 2024-12-03 [1] CRAN (R 4.3.0)
## rstudioapi 0.17.1 2024-10-22 [1] CRAN (R 4.3.0)
## S4Arrays 1.2.1 2024-03-04 [1] Bioconductor 3.18 (R 4.3.0)
## S4Vectors * 0.40.2 2023-11-23 [1] Bioconductor 3.18 (R 4.3.0)
## S7 0.2.0 2024-11-07 [1] CRAN (R 4.3.0)
## sass 0.4.9 2024-03-15 [1] CRAN (R 4.3.0)
## scales * 1.4.0 2025-04-24 [1] CRAN (R 4.3.0)
## scatterpie 0.2.4 2024-08-28 [1] CRAN (R 4.3.0)
## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.0)
## shadowtext 0.1.4 2024-07-18 [1] CRAN (R 4.3.0)
## shape 1.4.6.1 2024-02-23 [1] CRAN (R 4.3.0)
## SparseArray 1.2.4 2024-02-11 [1] Bioconductor 3.18 (R 4.3.0)
## statmod 1.5.0 2023-01-06 [1] CRAN (R 4.3.0)
## stringi 1.8.7 2025-03-27 [1] CRAN (R 4.3.0)
## stringr * 1.5.2 2025-09-08 [1] CRAN (R 4.3.0)
## SummarizedExperiment * 1.32.0 2023-10-24 [1] Bioconductor
## svglite 2.1.3 2023-12-08 [1] CRAN (R 4.3.0)
## systemfonts 1.2.1 2025-01-20 [1] CRAN (R 4.3.0)
## tibble * 3.3.0 2025-06-08 [1] CRAN (R 4.3.0)
## tidygraph 1.3.1 2024-01-30 [1] CRAN (R 4.3.0)
## tidyr * 1.3.1 2024-01-24 [1] CRAN (R 4.3.0)
## tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.3.0)
## tidytree 0.4.6 2023-12-12 [1] CRAN (R 4.3.0)
## tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.3.0)
## timechange 0.3.0 2024-01-18 [1] CRAN (R 4.3.0)
## treeio 1.26.0 2023-10-24 [1] Bioconductor
## tweenr 2.0.3 2024-02-26 [1] CRAN (R 4.3.0)
## tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.3.0)
## variancePartition * 1.32.5 2024-02-16 [1] Bioconductor 3.18 (R 4.3.0)
## vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.3.0)
## viridis 0.6.5 2024-01-29 [1] CRAN (R 4.3.0)
## viridisLite 0.4.2 2023-05-02 [1] CRAN (R 4.3.0)
## vroom 1.6.5 2023-12-05 [1] CRAN (R 4.3.0)
## withr 3.0.2 2024-10-28 [1] CRAN (R 4.3.0)
## xfun 0.53 2025-08-19 [1] CRAN (R 4.3.0)
## xml2 1.4.0 2025-08-20 [1] CRAN (R 4.3.0)
## XVector 0.42.0 2023-10-24 [1] Bioconductor
## yaml 2.3.10 2024-07-26 [1] CRAN (R 4.3.0)
## yulab.utils 0.1.9 2025-01-07 [1] CRAN (R 4.3.0)
## zlibbioc 1.48.2 2024-03-13 [1] Bioconductor 3.18 (R 4.3.0)
##
## [1] /opt/R/4.3.0/lib/R/library
##
## ──────────────────────────────────────────────────────────────────────────────